#include <graphics/perlinnoise.h>
#include <math.h>
#include <assert.h>


// Algorithms from:
// - http://freespace.virgin.net/hugo.elias/models/m_perlin.htm
// - http://www.sorgonet.com/linux/noise_textures/

namespace yam3d
{
namespace graphics
{

namespace
{
    static const float oneDiv4= 1.0f/4.0f;
    static const float oneDiv8= 1.0f/8.0f;
    static const float oneDiv16 = 1.0f/16.0f;
    static const float oneDiv32 = 1.0f/32.0f;
    static const float oneDiv64 = 1.0f/64.0f;
}


PerlinNoise::PerlinNoise( float persistence )
    : m_noiseFunc1(0)
    , m_noiseFunc2(0)
    , m_noiseFunc3(0)
    , m_interpolationFunc2(0)
    , m_interpolationFunc4(0)
    , m_persistence(0.0f)
    , m_width(128.0)
    , m_height(128.0)
{
    setPersistence(persistence);
    setSmoothing(false);
    setInterpolation(LINEAR);
}

PerlinNoise::~PerlinNoise()
{
}

void PerlinNoise::setPersistence(float persistence)
{
    m_persistence = persistence;
}

void PerlinNoise::setSmoothing(bool smoothed)
{
    if( smoothed )
    {
        m_noiseFunc1 = &PerlinNoise::noise1;
        m_noiseFunc2 = &PerlinNoise::noise2;
        m_noiseFunc3 = &PerlinNoise::noise3;
    }
    else
    {
        m_noiseFunc1 = &PerlinNoise::smoothedNoise1;
        m_noiseFunc2 = &PerlinNoise::smoothedNoise2;
        m_noiseFunc3 = &PerlinNoise::smoothedNoise3;
    }
}

void PerlinNoise::setInterpolation(InterpolationType interpolationType)
{
    m_interpolationFunc2 = 0;
    m_interpolationFunc4 = 0;
    switch(interpolationType)
    {
    case CLAMP:
        m_interpolationFunc2 = clampInterpolation;
        break;
    case LINEAR:
        m_interpolationFunc2 = linearInterpolation;
        break;
    case COSINE:
        m_interpolationFunc2 = cosineInterpolation;
        break;
    case CUBIC:
        m_interpolationFunc4 = cubicInterpolation;
        break;
    default:
        assert(0);
        break;
    }
}

float PerlinNoise::getNoiseValue(float x, int octave)
{
    if( m_interpolationFunc2 )
    {
        return perlinNoise(x, octave, m_interpolationFunc2 );
    }

    return perlinNoise(x, octave, m_interpolationFunc4 );
}

float PerlinNoise::getNoiseValue(float x, float y, int octave)
{
    if( m_interpolationFunc2 )
    {
        return perlinNoise(x, y, octave, m_interpolationFunc2 );
    }

    return perlinNoise(x, y, octave, m_interpolationFunc4 );
}

float PerlinNoise::getNoiseValue(float x, float y, float z, int octave)
{
    if( m_interpolationFunc2 )
    {
        return perlinNoise(x, y, z, octave, m_interpolationFunc2 );
    }

    return perlinNoise(x, y, z, octave, m_interpolationFunc4 );
}

float PerlinNoise::noise1(int32_t x)
{
  int32_t n = x;
  n = (n<<13) ^ n;
  return ( 1.0f - float( (n * (n * n * 15731 + 789221) + 1376312589) & 0x7fffffff) / 1073741824.0f);
}

float PerlinNoise::noise2(int32_t x, int32_t y)
{
  int32_t n = x + y*m_width;
  return noise1(n);
}

float PerlinNoise::noise3(int32_t x, int32_t y, int32_t z)
{
  int32_t n = x + y*m_width + z*m_width*m_height;
  return noise1(n);
}

float PerlinNoise::smoothedNoise1(int32_t x)
{
    return noise1(x)*0.5f  +  noise1(x-1)*0.25f  +  noise1(x+1)*0.25f;
}

float PerlinNoise::smoothedNoise2(int32_t x, int32_t y)
{
    // 4, weigth 1/16 0.25
    float corners = ( noise2(x-1, y-1) + noise2(x+1, y-1) + noise2(x-1, y+1) + noise2(x+1, y+1) );

    // 4, weigth 1/8 0.5
    float sides = ( noise2(x-1, y)   + noise2(x+1, y)   + noise2(x, y-1)   + noise2(x, y+1) );

    // 1, weigth 1/4 0.25
    float center = noise2(x, y);

    return oneDiv16*corners + oneDiv8*sides + oneDiv4*center;
}

float PerlinNoise::smoothedNoise3(int32_t x, int32_t y, int32_t z)
{
    float corners(0);
    float cornerSides(0);

    //8, weigth 1/64 = 0.125
    corners += ( noise3(x-1, y-1, z-1) + noise3(x+1, y-1, z-1) + noise3(x-1, y+1, z-1) + noise3(x+1, y+1, z-1) ); // z-1
    corners += ( noise3(x-1, y-1, z+1) + noise3(x+1, y-1, z+1) + noise3(x-1, y+1, z+1) + noise3(x+1, y+1, z+1) ); // z+1

    // 12, weigth 1/32 0.375
    cornerSides   += ( noise3(x-1, y, z)     + noise3(x+1, y, z)     + noise3(x, y-1, z)     + noise3(x, y+1, z) ); // z
    cornerSides   += ( noise3(x-1, y, z+1)   + noise3(x+1, y, z+1)   + noise3(x, y-1, z+1)   + noise3(x, y+1, z+1) ); // z+1
    cornerSides   += ( noise3(x-1, y, z-1)   + noise3(x+1, y, z-1)   + noise3(x, y-1, z-1)   + noise3(x, y+1, z-1) ); // z-1

    // 6, weigth 1/16 0.375
    float sides = ( noise3(x-1, y, z)   + noise3(x+1, y, z)   + noise3(x, y-1, z)   + noise3(x, y+1, z) + noise3(x, y, z-1)   + noise3(x, y, z+1) );

    // 1, weigth 1/8 0.125
    float center = noise3(x, y, z);

    return oneDiv64*corners + oneDiv32*cornerSides + oneDiv16*sides + oneDiv8*center;
}


float PerlinNoise::interpolatedNoise(float x, float (*interpolate)( float v0, float v1, float x ) )
{
     int32_t intX = int(x);
     float fractional = x - intX;


     float v0 = (this->*PerlinNoise::m_noiseFunc1)(intX + 0);
     float v1 = (this->*PerlinNoise::m_noiseFunc1)(intX + 1);
     //float v3 = (this->*noiseFunc)(intX + 2);

     return interpolate(v0, v1, fractional);
}

float PerlinNoise::interpolatedNoise(float x, float (*interpolate)( float v0, float v1, float v2, float v3, float x ) )
{
     int32_t intX = int(x);
     float fractional = x - intX;

     float v0 = (this->*PerlinNoise::m_noiseFunc1)(intX - 1);
     float v1 = (this->*PerlinNoise::m_noiseFunc1)(intX + 0);
     float v2 = (this->*PerlinNoise::m_noiseFunc1)(intX + 1);
     float v3 = (this->*PerlinNoise::m_noiseFunc1)(intX + 2);

     return interpolate(v0, v1 , v2, v3, fractional);
}

float PerlinNoise::interpolatedNoise( float x, float y, float (*interpolate)( float v0, float v1, float x ) )
{
    int32_t intX = int(x);
    float fractionalX = x - intX;

    int32_t intY = int(y);
    float fractionalY = y - intY;

    float vals[2] = {0.0f,0.0f};

    for( int i=0; i<2; ++i )
    {
        float vx0 = (this->*PerlinNoise::m_noiseFunc2)(intX + 0, y+i-1);
        float vx1 = (this->*PerlinNoise::m_noiseFunc2)(intX + 1, y+i-1);
        vals[i] = interpolate(vx0, vx1, fractionalX);
    }

    return interpolate(vals[0],vals[1],fractionalY);
}

float PerlinNoise::interpolatedNoise( float x, float y, float (*interpolate)( float v0, float v1, float v2, float v3, float x ) )
{
    int32_t intX = int(x);
    float fractionalX = x - intX;

    int32_t intY = int(y);
    float fractionalY = y - intY;

    float vals[4] = {0.0f,0.0f,0.0f,0.0f};

    for( int i=0; i<4; ++i )
    {
        float vx0 = (this->*PerlinNoise::m_noiseFunc2)(intX - 1, y+i-1);
        float vx1 = (this->*PerlinNoise::m_noiseFunc2)(intX + 0, y+i-1);
        float vx2 = (this->*PerlinNoise::m_noiseFunc2)(intX + 1, y+i-1);
        float vx3 = (this->*PerlinNoise::m_noiseFunc2)(intX + 2, y+i-1);
        vals[i] = interpolate(vx0, vx1 , vx2, vx3, fractionalX);
    }

    return interpolate(vals[0],vals[1],vals[2],vals[3],fractionalY);
}

float PerlinNoise::interpolatedNoise( float x, float y, float z, float (*interpolate)( float v0, float v1, float v2, float v3, float x ) )
{
    int32_t intX = int(x);
    float fractionalX = x - intX;

    int32_t intY = int(y);
    float fractionalY = y - intY;

    int32_t intZ = int(z);
    float fractionalZ = z - intZ;

    float valsZ[4] = {0.0f,0.0f,0.0f,0.0f};
    for( int i=0; i<4; ++i )
    {
        float valsY[4] = {0.0f,0.0f,0.0f,0.0f};
        for( int j=0; j<4; ++j )
        {
            float vx0 = (this->*PerlinNoise::m_noiseFunc3)(intX - 1, y+j-1, z+i-1);
            float vx1 = (this->*PerlinNoise::m_noiseFunc3)(intX + 0, y+j-1, z+i-1);
            float vx2 = (this->*PerlinNoise::m_noiseFunc3)(intX + 1, y+j-1, z+i-1);
            float vx3 = (this->*PerlinNoise::m_noiseFunc3)(intX + 2, y+j-1, z+i-1);
            valsY[j] = interpolate(vx0, vx1 , vx2, vx3, fractionalX);
        }

        valsZ[i] = interpolate(valsY[0],valsY[1],valsY[2],valsY[3],fractionalY);
    }

    return interpolate(valsZ[0],valsZ[1],valsZ[2],valsZ[3],fractionalZ);
}

float PerlinNoise::interpolatedNoise( float x, float y, float z, float (*interpolate)( float v0, float v1, float x ) )
{
    int32_t intX = int(x);
    float fractionalX = x - intX;

    int32_t intY = int(y);
    float fractionalY = y - intY;

    int32_t intZ = int(z);
    float fractionalZ = z - intZ;

    float valsZ[2] = {0.0f,0.0f};
    for( int i=0; i<2; ++i )
    {
        float valsY[4] = {0.0f,0.0f};
        for( int j=0; j<2; ++j )
        {
            float vx0 = (this->*PerlinNoise::m_noiseFunc3)(intX - 1, y+j-1, z+i-1);
            float vx1 = (this->*PerlinNoise::m_noiseFunc3)(intX + 0, y+j-1, z+i-1);
            valsY[j] = interpolate(vx0, vx1, fractionalX);
        }

        valsZ[i] = interpolate(valsY[0],valsY[1],fractionalY);
    }

    return interpolate(valsZ[0],valsZ[1],fractionalZ);
}

float PerlinNoise::powerOf2( int o )
{
    return pow(2,o);
}

float PerlinNoise::powerOfPersistence( int o )
{
    return pow(m_persistence,o);
}

float PerlinNoise::clampInterpolation( float a, float b, float x )
{
    return  x < 0.5f ? b : a;
}

float PerlinNoise::linearInterpolation( float a, float b, float x )
{
    return  a*(1.0f-x) + b*x;
}

float PerlinNoise::cosineInterpolation( float a, float b, float x )
{
    float ft = x * 3.1415927f;
    float f = (1.0f - cosf(ft)) * 0.5f;
    return  a*(1.0f-f) + b*f;
}

float PerlinNoise::cubicInterpolation( float v0, float v1, float v2, float v3, float x )
{
    float P = (v3 - v2) - (v0 - v1);
    float Q = (v0 - v1) - P;
    float R = v2 - v0;
    float S = v1;
    return P*(x*x*x) + Q*(x*x) + R*x + S;
}



}
}

